\documentclass[9pt,twocolumn,twoside]{osajnl}
%% Please use 11pt if submitting to AOP
% \documentclass[11pt,twocolumn,twoside]{osajnl}

\journal{optica} % Choose journal (ao, aop, josaa, josab, ol, optica, pr)

\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{amsmath}
\usepackage{cleveref}

\DeclareMathOperator{\sinc}{sinc}
\DeclareMathOperator{\somb}{somb}
% See template introduction for guidance on setting shortarticle option
\setboolean{shortarticle}{false}
% true = letter / tutorial
% false = research / review article
% (depending on journal).
\usepackage[justification=centering,
            format=plain]{caption}

\title{Paper}

% \author[1,2]{Haotian Song}
% \author[1,3,4,$\dagger$]{Xiaoyu Nie}
% \author[1]{Zhaoyu Zuo}
% \affil[1]{Physics Department, Xi'an Jiaotong University, Xi'an, Shaanxi 710049, China}
% \affil[2]{School of Physics \& Astronomy, University of Manchester, Manchester M13 9PL, United Kingdom}
% \affil[$\dagger$]{nxy1125@tamu.edu}
\setboolean{displaycopyright}{true}
\begin{abstract}
This PDF show the process of caculating SFR and $\delta N$
 
\end{abstract}
\begin{document}
\maketitle
\section{INTRODUCTION}



\section{Methods}

We utilized the Binary Population Synthesis (BPS) code \cite{hurley2002evolution} with further updates, which are briefly described as follows. And MESA was performed in the population synthesis of NS for more specific evolution. 

% We set-up a grid of initial binary parameters as following.
% \begin{equation}\begin{aligned}
%     M_{1} &: \quad 2 \rightarrow 150 \mathrm{M}_{\odot} \\
%     q &: \quad 0.08 \rightarrow 1 \mathrm{M}_{\odot} \\
%     a &: \quad 3.0 \rightarrow 10^{4} \mathrm{R} \odot
%     \end{aligned}
% \end{equation}
The Initial Mass Function (IMF) developed by Kroupa, Tout $\&$ Gilmore (1993) \cite{kroupa1993distribution} was applied. Although the shapes of IMF differ significantly, it is less subject to the uncertainty due to the similarity of the IMFs for larger masses. The mass of star on zero-age sequence (ZAMS) satisfies the following distribution with initial mass ranging from 2 $M_\odot$ to 150 $M_\odot$.
\begin{equation}
    \xi(m) \propto m^{-\alpha}
\end{equation}
where
\begin{equation}
    \alpha(m)=\left\{\begin{array}{ll}
        +0.3 \pm 0.7 & 0.01 \leq m < 0.08 \\
        +1.3 \pm 0.5 & 0.08 \leq m < 0.50 \\
        +2.3 \pm 0.3 & 0.50 \leq m < 1.00 \\
        +2.3 \pm 0.7 & 1.00 \leq m
           \end{array}\right.
\end{equation}
The distribution of mass ratios q ($q=M_1/M_2$) was uniform between 0.08 and 1 by steps of 0.10 \cite{wiktorowicz2019observed}. The orbital separation was ranged from 3.0 to $10^4 R_\odot$ which facilitates the comparison with other models \cite{yungelson1997rate}. Here we assumed that $ln(a)$ is uniformly distributed. The maximum evolution time was 200Myr in isolation for comparison with young galaxies. As for supernova kicks, we draw them from a Maxwellian distribution with $\sigma=265km \cdot s^{-1}$ \cite{hobbs2005statistical}. And solar metallicity $Z=0.5Z_\odot$ was adopted in our simulations. 

For sources undergoing Roche Lobe Overflow (RL) mass-transfer (MT), the traditional formula was used for sub-Eddington accretion rates, where Eddington luminosity is $L_E = 1.6 \times 10^{38} m_1$ for hydrogen-rich material. Supplied with super-Eddington mass transfer rates, the accretor expel matter in significantly different approach. Outside the spherization $R_{sph}$ \cite{blundell2007fluctuations}, accretion luminosity is released as usual. \cite{shakura1973black} \cite{poutanen2007supercritically} But within $R_{sph}$, radiation become inefficient and the outflow maintain energy release in $L_E$ under the pressure of radiation. Bolometric luminosity $L_X$ is then
\begin{equation}L_{X}=\left\{\begin{array}{ll}
    L_{\mathrm{Edd}}\left(1+\ln \dot{m}_{\mathrm{tr}}\right) & \dot{m}_{\mathrm{tr}}>1 \\
    L_{\mathrm{Edd}} \dot{m}_{\mathrm{tr}} & \dot{m}_{\mathrm{tr}} \leqslant 1
    \end{array}\right.
\end{equation}
Where $\dot{m}_{\mathrm{tr}}=\dot{M}_{\mathrm{tr}} / \dot{M}_{\mathrm{Edd}}$ and $\dot{M}_{\mathrm{tr}}$ is in units of $M_\odot/yr$. For the region within $R_{sph}$, a biconical geometry is formed with collimation of radiation\cite{king2008accretion}. Therefore, the finite cones can be detected by observer. The apparent (isotropic) X-ray luminosity is 
\begin{equation}
    L_{app}=L_X/b
\end{equation}

Where b is beaming factor defined as the fraction of funnels area, $b \stackrel{\text { def }}{=} \Omega / 4 \pi$, and $\Omega$ is the combined solid angle of both beams. When $\dot{m}$ is larger enough, the scales can vary from the typical cylinderical radius to $R_{sph}$  markedly. And for soft-excess luminosity, the observed relation $L_{\mathrm{soft}} \propto T^{-4}$ sustain the estimated value of b \cite{king2009masses}.
\begin{equation}
    b=\left\{\begin{array}{ll}
    1 & \dot{m}_{\mathrm{tr}} \leqslant 8.5 \\
    \frac{73}{\dot{m}_{\mathrm{tr}}^{2}} & 8.5<\dot{m}_{\mathrm{tr}}
    \end{array}\right.
\end{equation}
whereas for $ \dot{m} < 8.5$, it is assumed unbeamed. And this formula have been verified in plenty of hyperluminous X-ray sources \cite{king2014hlx,king2016ulxs,king2017pulsing}. In addition, there was an implication \cite{wiktorowicz2017origin} that extremely beamed sources are scarced and transient that lead to difficulty in detection. That means there is no lower valve for the beaming factor.

Furthermore, the wind Roche lobe overflow (WRL) could be taken into acount. The traditional Bondi-Hoyle-Lyttleton (BHL) MT formula for WRL \cite{bondi1944mechanism,edgar2004review} is 

\begin{equation}
    \dot{M}_{\mathrm{BHL}}=\pi R_{\mathrm{BHL}}^{2} v_{\mathrm{rel}} \rho
\end{equation}

where $\dot{M}_{\mathrm{BHL}}$ is the mass accretion rate, $v_{\mathrm{rel}}=\sqrt{v_{\beta}^{2}+(v_{\text {orb}}[q /(1+q)])^{2}}$ is relative speed between the wind and compact object. $v_{\text {orb}}$ is orbital speed given by $v_{orb}=2\pi a / P_{orb}$ with $P_{orb}$ being the orbit period. $ R_{BHL}=2 G M_{\bullet} / v_{\mathrm{rel}}^{2}$ is the modified accretion radius and $\rho$ is the density of the wind at orbital separation. After assumption of isotropic dilution of the wind, the fraction of wind captured $\mu_{\mathrm{BHL}}=\dot{M}_{\mathrm{BHL}} / \dot{M}_{\star}$ can be:

\begin{equation}
    \mu_{\mathrm{BHL}}=\frac{(1+q) / q^{3}}{\eta(1-f \mathcal{E})^{\beta}\left[1+\left(\eta(1+q)(1-f \mathcal{E})^{\beta} / q\right)\right]^{3 / 2}}
\end{equation}

where $\mathcal{E}$ is the ratio of stellar Roche lobe radius by the orbital separation and only related to q \cite{eggleton1983aproximations}, and $\eta$ is the speed ratio. However, another research \cite{el2019wind} proposed that ULXs could remain stable in WRL with a highly beamed wind. Based on realistic acceleration profiles, they estimated the MT by computing the bulk motion of the wind. By assuming non-isotropic dilution of stellar wind, the fraction of wind captured can be significantly higher especially for higher wind speed ratio. And several binary systems have been proved possible for WRL, such as M101 \cite{liu2013puzzling}, and P13 \cite{furst2018tale}. We adopted his simulation and assume $\mu$ as a function of $\eta$, stellar filling factor, $\beta$ and mass ratio $q$. The WRL efficiency systematically improved. And WRL can be a possible mechanism for ULXs.

As for the observation of ULXs, the statistics of ULXs in seven galaxies have been published \cite{wolter2018x}. 
%Relevant properties are listed in Table???. 
For numerical comparison with them, the evolved population of binaries should be in conjunction with a realistic birth rates and metallicity. Here the distribution of each binary system is \cite{hurley2002evolution}
 \begin{equation}
    \label[]{eq:hurley}
    \delta r_{j}=\overline{S_b} \Phi\left(\ln M_{1 j}\right) \varphi\left(\ln M_{2 j}\right) \Psi\left(\ln a_{j}\right) \delta \ln M_{1} \delta \ln M_{2} \delta \ln a
\end{equation}
where r is the rate of particular population, and $\overline{S_b}$ is the binary star forming rate (SFR) in units of $number/year$ which can be obtained from following formula with SFR in units of $M_\odot/yr$ for binaries massive than 5 $M_\odot$ \cite{grimm2003high}:
\begin{equation}
    \overline{S_b}=(\frac{M_{low}}{5})^{1-\alpha}\frac{(\alpha-2)SFR}{5(\overline{3}+\overline{q})(\alpha-1)}
\end{equation}
where $M_{low}$ is the lower limit initial mass for BPS. For summation of seven galaxies $SFR\approx 42 M_\odot/yr$, $\overline{S_b}=1.5039 yr^{-1}$ The total number N should be the summation of every binary system distribution. However, the number observed could differ in consideration of beaming effect. With assumation of linear relation between the area of the flux sphere and beaming factor, the chance probability is therefore $P \propto b$ \cite{Middleton2017}. The observed number $N_{OBS}$ is 
\begin{eqnarray}
    N_{OBS}=\sum \delta N = b \sum \delta r \delta T 
\end{eqnarray}
where $ \delta T $ is the duration of particular binary evolution stage. For mesa simulation, $\delta r$ was obtained from BPS code where particular event is the birth of NS.

\section{Simulation Results}

The simulation in this section was performed with constant SFR and metallicity. Maxium evolution time is 200Myr for early formed starburst galaxies. And the number of ULXs is linear to constant SFR in this algorithm which accords with previous work \cite{Mineo2013}. With assumption of constant SFR for observed galaxies, we compare it with observation and previous analysis. Due to complex episodes of star forming, estimates should be modified via a more realistic method. 

Starburst galaxies are comparatively rare (15 out of 345 galaxies), in which, however, the ULXs occupy the certain proportion \cite{Wang2016}. As a concequece, starburst galaxies can imply general characteristic of ULXs. We have conducted our simulation to compare with galaxy-focused observations \cite{wolter2018x}. The result of luminosity function is presented in Fig. \ref{fig:mainLF}.

\begin{figure}[ht!]
    % \centering
    \subfigure[]{
		\begin{minipage}[t]{\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/MAINLuminosityFunction.jpg}
            
		\end{minipage}
    }%
    
    \subfigure[]{
		\begin{minipage}[t]{\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/LuminosityFunctionH-He.jpg}
		\end{minipage}
	}%
	\caption{X-ray luminosity function(XLF) for seven starburst galaxies mentioned above. The solid blue line represents total ULXs number via BPS simulation, and red dots represent observation. (a) Dashed dotted lines represent the ULXs number in Roche Lobe accretion and wind accretion respectively. The compact star is Black Hole(BH) and Neutron Star(NS) for pale blue and green lines respectively. (b) The simulation distribution of ULX donor and compact star. Dashed and dotted lines are ULX with H and He donor star respectively, and orange lines are BH compact star while green line is NS.}
	\label{fig:mainLF}
\end{figure}

The estimated numbers of ULXs are based on Eq. \ref{eq:hurley}. The summation of the SFR for 7 galaxies is approximately $43.5M_\odot$ obtained from $L_{H\alpha}$ \cite{appleton1997multiwavelength}. Large uncertainty of measuring SFR should be noticed. Therefore the observation of limited galaxies may not be representive and comparable to our simulation. 50 sources above the ULX valve ($10^{39}erg\cdot s^{-1}$ \cite{Kaaret2017}) are collected of which 23 sources are above $5 \times 10^{39}erg\cdot s^{-1}$. The BPS simulation generally fits the observation and display the characteristic of ULXs. The pure power-law sample was presented for XLF \cite{Swartz2011} where 107 ULXs were identified. But the census of these galaxies is not suitable for our simulation for starburst galaxies. 
%  The relation can be expressed as $C L_{\mathrm{X}}^{-\alpha} \exp \left(-L_{\mathrm{X}} / L_{c}\right)$ where $L_C$ is the cutoff luminosity, C is the normalization constant, and $\alpha$ is the power-law index. With the minimum-likelihood fit statistic, the parameters are $C = 103.3, \alpha = 0.9$ in his work. 

According to the BPS simulation result, there is a break at approximately $3 \times 10^{40} erg \cdot s^{-1}$ which accords with the observation of other galaxies \cite{mineo2012x,Swartz2011}. Although extreme objects reaching a peak luminosity of $10^{41} erg \cdot s^{-1}$ exist (e.g., M82) \cite{gao2003nonnuclear}, they are candidates for Intermediate Mass Black Holes (IMBHs) with masses in between ($ 10^2 ~ 10^5 M_\odot $) \cite{kuranov2007dynamical}. With completely disparate formation scenarios, IMBHs exclude from consideration. Extremely luminous sources ($L_X > 3 \times 10^{40} erg \cdot s^{-1} $) demand high beaming factor and correspondingly are impossible to be detected. 
For less luminous sources, the compact stars of ULXs tend to be BH instead of NS. The typical orbit period for BH via RL is a few days and donor star mass tends to be less than 2.7 solar mass. Meanwhile, NS dominates in the region with $L_{app} > 10 ^ {40} erg \cdot s^{-1}$, which is analogous to the previous work \cite{Shao2015}. The accretion via wind Roche lobe overflow contributes quiet a few by implying new mechenism.

\begin{figure}[ht!]
    \includegraphics[width=\linewidth]{pictures/TimeTotal.jpg}
    \caption{The evolution duration of ULXs for 200 Myr, 100 Myr, 50 Myr, 20 Myr, 15 Myr.}
    \label{fig:TimeXLF}
\end{figure}

In Fig. \ref{fig:TimeXLF}, various evolution duration was performed in the simulation. The number of more luminous ULXs significantly drops at an age of ~150Myr whose compact star corresponds to NS. When the population is around 30 Myr, ULXs harboring BH dominate. Galaxies younger than 20 Myr contributes a few because this requires strict condition for formation of compact star.

\begin{figure}[ht!]
    \includegraphics[width=\linewidth]{pictures/Lx-Porb-total.jpg}
    \caption{The numerical distribution in the $P_{orb}-L_X$ plane for ULXs. The depth of shade represents the quantity.}
    \label{fig:LxPorb-total}
\end{figure}

Fig. \ref{fig:LxPorb-total} presents the numerical distribution of ULX with luminosity larger than $10^{39}erg\cdot s^{-1}$. It is seen that binaries with orbit periods shorter than few days dominate, which corresponds with previous work \cite{Shao2015}. Fig. \ref{fig:M2Porb-total} is the same with Fig. \ref{fig:LxPorb-total} but for $M_2-L_X$ plane. Due to larger possibility in IMF, lower massive donor star comprise a large proportion.
% In Fig. \ref{fig:LxPorbIndividual} and Fig. \ref{fig:M2PorbIndividual}, various types of ULXs are displayed individually for their orbit periods, luminosity and masses of donor stars.

\begin{figure}[ht!]
    \includegraphics[width=\linewidth]{pictures/M2-Porb-total.jpg}
    \caption{Same as Fig. \ref{fig:LxPorb-total}, but for $M_2-L_X$ plane.}
    \label{fig:M2Porb-total}
\end{figure}

% \begin{figure*}[ht!]
	
% 	\subfigure{
% 		\begin{minipage}[t]{0.5\linewidth}
% 			\centering
% 			\includegraphics[width=\linewidth]{Lx-Porb-BH-RLOF.jpg}\includegraphics[width=\linewidth]{Lx-Porb-BH-wind.jpg}
% 			\includegraphics[width=\linewidth]{Lx-Porb-NS-RLOF.jpg}\includegraphics[width=\linewidth]{Lx-Porb-NS-wind.jpg}
% 			% \caption{BH-Roche Lobe}
% 		\end{minipage}
% 	}
% 	\centering
% 	\caption{Predicted census of ULXs for several specific types. In the upper and lower panels, the compact stars are BH and NS respectively. In the left and right panels, the accretion mechenism is Roche Lobe and wind respectively.}
% 	\label{fig:LxPorbIndividual}
% \end{figure*}


% \begin{figure*}[ht!]
	
% 	\subfigure{
% 		\begin{minipage}[t]{0.5\linewidth}
% 			\centering
% 			\includegraphics[width=\linewidth]{M2-Porb/M2-Porb-BH-RLOF.jpg}\includegraphics[width=\linewidth]{M2-Porb/M2-Porb-BH-wind.jpg}
% 			\includegraphics[width=\linewidth]{M2-Porb/M2-Porb-NS-RLOF.jpg}\includegraphics[width=\linewidth]{M2-Porb/M2-Porb-NS-wind.jpg}
% 			% \caption{BH-Roche Lobe}
% 		\end{minipage}
% 	}
% 	\centering
% 	\caption{Same as Fig. \ref{fig:LxPorbIndividual}, but for $P_{orb}-M2$ plane.}
% 	\label{fig:M2PorbIndividual}
% \end{figure*}






%
	% \subfigure[]{
	% 	\begin{minipage}[t]{0.5\linewidth}
	% 		\centering
	% 		\includegraphics[width=\linewidth]{M2-Porb/M2-Porb-BH-wind.jpg}
	% 		% \caption{BH-Wind}
	% 	\end{minipage}
    % }%
    
    %此处的空行很重要，想让图片在什么地方换行就在代码对应位置空行
	% \subfigure[]{
	% 	\begin{minipage}[t]{0.5\linewidth}
	% 		\centering
	% 		\includegraphics[width=\linewidth]{M2-Porb/M2-Porb-NS-RLOF.jpg}\includegraphics[width=\linewidth]{M2-Porb/M2-Porb-NS-wind.jpg}
	% 		% \caption{NS-Roche Lobe}
	% 	\end{minipage}
	% }%
	% \subfigure[]{
	% 	\begin{minipage}[t]{0.5\linewidth}
	% 		\centering
	% 		\includegraphics[width=\linewidth]{M2-Porb/M2-Porb-NS-wind.jpg}
	% 		% \caption{NS-Wind}
	% 	\end{minipage}
	% }%
% The following is specific caculation for 
% \begin{equation}
%     L_{\mathrm{X}}=\left\{\begin{array}{ll}
%     L_{\mathrm{Edd}}\left(1+\ln \dot{m}_{\mathrm{tr}}\right) & \dot{m}_{\mathrm{tr}}>1 \\
%     L_{\mathrm{Edd}} \dot{m}_{\mathrm{tr}} & \dot{m}_{\mathrm{tr}} \leqslant 1
%     \end{array}\right.
% \end{equation}

% where $\dot{m}_{tr}=\dot{M}_{tr}/\dot{M}_{Edd}$
% \begin{equation}
%     L_{edd}=\eta \dot{M}_{Edd} c^2
% \end{equation}
% whether $\eta$ is 1 or 0.1,


% I am not sure whether SFR includes all stars or binaries only, so I caculate it individually.

% $S_b$and $S_s$ are in units of Num/yr.

% SFR is in units of $M_\odot/yr$

% The main idea of caculation is as following.
% \begin{equation}
%     \label{eq:1}
%     S_b \dot (average of M_{low} and M_{up}) + S_b \dot (average M_2) + S_s \dot(average mass) = SFR
% \end{equation}

% Equation\ref{eq:1}can be written as following, too.

% \begin{equation}
%     S_b \dot \int_{M_{low}}^{M_{up}} m_1 \epsilon({m_1})dm_1 + S_b \int_{M_{low}}^{M_{up}} \epsilon({m_1}) (\int_{M_{low}}^{M_1}\frac{m_2}{m_1}  dm_2) dm_1 + S_s \int_{M_{low}}^{M_{up}} m \epsilon({m})dm = SFR
% \end{equation}


% \begin{equation}
%     S_b \dot \int_{M_{low}}^{M_{up}} m_1 \epsilon({m_1})dm_1 + S_b \int_{M_{low}}^{M_{up}} \epsilon({m_1}) m_1 \overline{q}  dm_1 + S_s \int_{M_{low}}^{M_{up}} m \epsilon({m})dm = SFR
% \end{equation}

% chosen q in my caculation is 0.08, 0.18, 0.28... 0.98, so $\overline{q}=0.535$

% If I chose $M_{low}=5M_\odot$, as\cite{grimm2003high} .
% \begin{equation}
%     SFR=(S_b(1+\overline{q}) + S_s)\int_{M_{low}}^{\infty} m \epsilon(m)dm
% \end{equation}

% However, the Initial Mass Function differs.
% When $m>1M_\odot$, $\epsilon(m)=am^{-\alpha}$

% For \cite{kroupa2001variation}, $\alpha=+2.3\pm0.7$

% For \cite{kroupa1993distribution} $\alpha=+2.7$

% \begin{equation}
%     SFR=(S_b (1+\overline{q}) + S_s)\int_{M_{low}}^{\infty} am^{1-\alpha}dm
% \end{equation}
% \begin{equation}
%     SFR=(S_b (1+\overline{q}) + S_s) a\frac{m^{2-\alpha}}{2-\alpha}\bigg|_{M_{low}}^{\infty}
% \end{equation}
% \begin{equation}
%     SFR=(S_b (1+\overline{q}) + S_s) a\frac{M_{low}^{2-\alpha}}{\alpha-2}
% \end{equation}

% For my BPS caculation, $M_{low}=2M_\odot$, so
% \begin{equation}
% \overline{S_b}=S_b\frac{\int_2^\infty\epsilon(m)dm}{\int_5^\infty\epsilon(m)dm}=S_b(\frac{2}{5})^{1-\alpha}
% \end{equation}

% And the relation between $S_b$ and $S_s$ is 
% $\frac{2S_b}{S_s+2S_b}=f=0.5$
% where f is the fraction of binaries and all stars. If $f=0.5$, then $S_s=2S_b$

% And we will get a in the following.

% As is shown before,  $\epsilon(m)=am^{-\alpha}$. And  $\epsilon(m)$ should be normalized from $M_{low}$ to $M_{max}$, so a is changed and not equal to the Initial a in the essay.
% $$1=\int_{M_{low}}^{M_{max}}am^{-\alpha}dm$$

% $$1=a\frac{m^{1-\alpha}}{1-\alpha}\bigg|_{M_{low}}^{\infty}=a\frac{M_{low}^{1-\alpha}}{1-\alpha}$$

% $$a=\frac{\alpha-1}{M_{low}^{1-\alpha}}$$

% In conclution, $S_b$ can be written as following.
% \begin{equation}
%     \overline{S_b}=(\frac{2}{5})^{1-\alpha}\frac{(\alpha-2)SFR}{(\overline{3}+\overline{q})M_{low}(\alpha-1)}
% \end{equation}
% where $\overline{3}$ means whether SFR contains single star. If not, $\overline{3}=1$.

% \begin{table}[h!]
% \centering
% \begin{tabular}{lll}
%          & $\alpha=2.3$ & $\alpha=2.7$ \\
% Binary   & 4.1560       & 10.6984      \\
% All star & 1.8046       & 4.6455
      
% \end{tabular}
% \end{table}
% %\noindent\textbf{Disclosures.} The authors declare no conflicts of interest.
% \section{Caculate $\delta N$ }
% According to \cite{hurley2002evolution}, the $\delta r$ can be written as:
% \begin{equation}
% \delta r_{j}=S \Phi\left(\ln M_{1 j}\right) \varphi\left(\ln M_{2 j}\right) \Psi\left(\ln a_{j}\right) \delta \ln M_{1} \delta \ln M_{2} \delta \ln a
% \end{equation}
% and $\delta N = \delta r \delta T$, so
% \begin{equation}
%     \delta N = \overline{S_b} M_1 \epsilon({m_1}) \delta \ln M_1 k_q \delta q \dot k_a \delta \ln a \delta T
% \end{equation}
% And the q and ln(a) are mean distribution, so $k_q\delta q = 1/(q_{number})$,$k_a\delta \ln a=1/a_{number}$
% where $q_{number}$ and $a_{number}$ is the cycle index in my caculation. $1/a_{number}$ is the length of step after normalization.
% \begin{equation}
%     \delta N = \overline{S_b} a M_1^{1-\alpha} \delta \ln M_1 k_q \delta q \dot k_a \delta \ln a \delta T
% \end{equation}

% \begin{equation}
%     \delta N = \overline{S_b} (\alpha-1) (\frac{M_1}{M_low})^{1-\alpha} \frac{\ln(M_{max})-\ln(M_{min})}{M_{number}} \frac{1}{q_{number}} \frac{1}{a_{number}} \delta T
% \end{equation}










% \section{Adding A to the formula}

% The sum of all Success cases should be 1, 
% \begin{equation}
%     \Sigma Success  = 1,\iiint M_1 \epsilon({m_1}) \delta \ln M_1 k_q \delta q \dot k_a \delta \ln a=1
% \end{equation}
% but it is not because of failure cases. So A is added.
% \begin{equation}
%     \Sigma Success + \Sigma Failure  = A \iiint M_1 \epsilon({m_1}) \delta \ln M_1 k_q \delta q \dot k_a \delta \ln a 
% \end{equation}
% \begin{equation}
%     \Sigma Success + \Sigma Failure  = 1 + \Sigma Failure
% \end{equation}
% \begin{eqnarray}
%     A=\frac{1}{ \iiint - \Sigma Failure}=A
% \end{eqnarray}

% \begin{equation}
%     \delta N =(1+\Sigma Failure) \overline{S_b} (\alpha-1) (\frac{M_1}{2})^{1-\alpha} \frac{\ln(M_{max})-\ln(M_{min})}{M_{number}} \frac{1}{q_{number}} \frac{1}{a_{number}} \delta T
% \end{equation}

% And I got $A=1.0295$
\bibliography{sample}


\end{document}